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INTRODUCTION 


In  the  high  temperature  environment  of  advanced  high  performance 
turbines  of  interest  to  the  Air  Force,  cooling  air  and  its  associated  cost 
becomes  a  critical  item.  Typically  cooling  air  is  bled  from  the  compressor 
and  delivered  to  the  turbine  through  the  machine  core  out  of  the  primary  gas 
path  and  circumventing  the  combustor.  Provided  the  cooling  air  is  discharged 
into  the  turbine  sufficiently  early,  the  major  losses  are  incurred  in  the  off 
gas  path  flow  and  delivery  system.  These  losses  can  be  an  appreciable  cost 
in  overall  efficiency  when  the  cooling  mass  flow  requirements  become  a 
significant  fraction  of  the  primary  gas  flow,  as  occurs  in  most  high 
temperature  high  performance  engines.  Taking  the  bleed  flow  out  of  the  gas 
path  and  delivering  it  to  the  turbine  with  little  loss  in  total  pressure 
turns  out  to  be  a  difficult  task.  Furthermore,  the  cooling  air  can  be 
required  also  to  cool  the  turbine  disks  prior  to  or  in  addition  to,  cooling 
the  turbine  rotor  or  vane  blades. 

A  major  problem  that  can  arise  in  this  cooling  air  delivery  process  is 
associated  with  ability  of  the  rotor  disk  to  pump  the  cooling  air  out  into 
the  gas  path.  The  so-called  'disk  pumping  problem'  arises  in  part  by  virtue 
of  the  shear  stress  on  the  rotating  disk.  This  is  in  itself  a  drag  producing 
mechanism.  Enough  air  is  required  to  be  driven  out  in  the  gas  path  to 
prevent  ingestion  of  the  hot  main  gas  flow  into  the  disk  cavity.  However, 
the  amount  involved  should  not  be  so  much  as  to  adversely  thicken  the  endwall 
boundary  layers  or  to  otherwise  disrupt  the  primary  gas  flow.  Furthermore, 
it  is  desirable  on  an  overall  basis  to  minimize  the  coolant  mass  flow 
requirements . 

Further  complications  can  arise  by  virtue  of  modern  turbine  design  which 
take  advantage  of  the  structural  savings  which  result  from  hanging  the 
stators  from  the  casing,  as  is  shown  schematically  in  Fig.  1.  The  axial 
static  pressure  drop  in  the  turbine  requires  the  stator  to  be  sealed  and  one 
possible  method  is  shown  in  Fig.  1.  Under  certain  conditions  the  resulting 
configuration  may  allow  ingestion  of  the  hot  main  flow  gas  into  the  cavity 
such  as  in  cavity  "d",  and  special  care  is  required  to  inhibit  or  prevent 
this  process.  Clearly  the  problems  of  disk  pumping,  disc  cooling  and  overall 
interaction  with  the  main  gas  path  are  complex.  In  view  of  their  importance 
to  structural  integrity,  durability  and  efficiency  it  would  be  very  desirable 
to  have  a  good  understanding  and  predictive  capability  for  these  flows. 


Insofar  as  Che  analysis  is  concerned,  unfortunately  the  complex 
turbulent  recirculating  rotational  nature  of  the  flow  does  not  permit 
significant  simplification  in  developing  a  mathematical  approach,  and 
consequently  to  date  very  little  has  been  done  in  the  way  of  developing 
formal  analyses  for  this  problem  from  the  basic  governing  equations.  Various 
geometry  dependent  variations  of  the  disk  pump  concept  can  arise,  however, 
only  a  limited  amount  of  experimental  data  for  these  various  possible 
configurations  is  available.  Hence,  an  analysis  capable  of  accurately  and 
efficiently  predicting  disk  cavity  flows  would  be  a  major  tool  for  both  the 
research  and  design  engineers.  The  initial  development  of  such  a  program  was 
the  major  focus  of  the  present  Phase  I  effort. 

Considering  that  the  disk  pump  may  operate  under  the  possible 
recirculating  and  turbulent  flow  conditions,  it  is  clear  that  the  solution  of 
the  ensemble-averaged  Navier-Stokes  equations  represent  the  appropriate 
approach  to  the  problem.  However,  accurate  and  efficient  solutions  of  the 
Navier-Stokes  equations  in  complex  geometries  for  realistic  flow  conditions 
represent  a  formidable  task.  Recently,  Buggeln  and  McDonald  (Ref.  1)  have 
successfully  solved  the  Navier-Stokes  equations  for  a  variety  of  labyrinth 
seal  conf igurat ions  at  realistic  flow  conditions,  and  have  obtained  results 
which  show  good  agreement  with  experimental  measurements.  Although  the  disk 
pumping  problem  differs  from  the  labyrinth  seal  problem,  it  was  felt  that 
this  same  procedure  with  modification  could  be  used  as  an  accurate  and 
efficient  approach  to  the  disk  pumping  problem.  Therefore,  under  the  present 
Phase  I  effort,  the  work  focused  upon  demonstrating  and  assessing  the 
applicability  of  a  modified  version  of  this  labyrinth  seal  computer  code  to 
the  problem  of  the  flow  in  a  disk  pumping  cavity  of  type  d  in  Fig.  1. 

In  the  present  study  an  axisymmetric  configuration  was  supposed  for  the 
disk  pump,  and  this  is  considered  a  good  approximation  for  the  present 
program  since  the  axisymmetric  configuration  gives  rise  to  the  relevant 
physical  flow  phenomena.  It  should  be  noted  that  the  axisymmetric 
configuration  permits  swirling  cavity  flows  to  arise,  but  would  not  treat  a 
discretely  three-dimensional  proturberance ,  such  as  a  bolt  head.  Fully 
three-dimensional  flows,  although  possible  to  compute  with  this  type  of 
analysis  (e.g.  Refs.  2  and  3),  would  require  considerably  more  computer  run 
time  than  axisymmetric  flows.  Therefore,  since  the  present  effort  was  to 
concentrate  upon  assessing  the  potential  of  the  approach,  it  was  decided  not 


to  pursue  three-dimensional  flows  under  Phase  I  of  the  effort. 
Three-dimensional  flows  as  well  as  other  complicating  factors  could  be 
performed  on  a  Phase  II  effort. 

In  order  to  demonstrate  the  ability  to  analyze  simple  disk  pumping 
flows,  the  present  effort  concentrated  on  predicting  the  flow  in  a  cavity  of 
the  type  d  shown  in  Fig.  1.  In  this  type  of  cavity  one  disk  wall  is  fixed, 
the  other  rotating.  Air  enters  along  the  disk  inner  radius  wall  and  leaves 
along  the  disk  periphery.  The  final  model  disk/cavity  configurations 
actually  utilized  in  the  present  study  are  shown  in  Fig.  2.  Although 
somewhat  simplified,  it  represents  a  viable  approximation  to  actual 
geometries . 

This  report  presents  and  describes  the  analysis  of  disk  pumping  flows 
based  on  the  linearized  block  implicit  (LBI)  technique  for  the  numerical 
solution  of  the  multidimensional,  time-dependent  Navier-Stokes  equations  by 
Briley  and  McDonald  (Ref.  4).  Two  different  shapes  of  disk  pump  were 
considered  in  the  demonstration  calculation  (Fig.  2).  In  the  following 
sections,  a  general  description  of  the  present  analysis,  procedures,  results 
and  an  assessment  of  the  potential  of  the  method  are  given. 


ANALYSIS 

Governing  Equations 


The  equations  used  in  the  present  effort  are  the  ensemble-averaged, 
time-dependent  Navier-Stokes  equations  which  can  be  written  in  vector  form  as 
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where  p  is  density,  u  is  velocity,  p  is  pressure,  it  is  the  molecular  stress 
tensor,  it1  is  the  turbulent  stress  tensor,  h  is  enthalpy,  q  is  the  mean 
heat  flux  vector,  q^  is  the  turbulent  heat  flux  vector,  ♦  is  the  mean  flow 
dissipation  rate  and  e  is  the  turbulence  energy  dissipation  rate.  If  the 
flow  is  assumed  to  have  a  constant  total  temperature,  the  energy  equation  is 
replaced  by 

2 

T  =  T  +  =  constant  (4) 

o  ZL 

P 

where  T0  is  the  stagnation  temperature,  q  is  the  magnitude  of  the  velocity 
and  Cp  is  the  specific  heat  at  constant  pressure.  In  both  cases  considered 
in  this  work  the  assumption  of  constant  total  temperature  has  been  invoked  by 
using  Eq.  4  as  an  approximation  to  Eq.  3  for  the  sole  purpose  of  reducing 
computer  run  time  where  the  constant  T0  assumption  was  warranted. 

A  number  of  terms  appearing  in  Eqs .  1-4  require  definition.  The  stress 
tensor  appearing  in  Eq.  2  is  defined  as 

*  -  2u  -  li-Kjj )  7  •  U  n  (5) 

where  Kg  is  the  bulk  viscosity  coefficient  and  is  the  deformation  tensor, 
defined  by: 

CD  =  ((VU)  +  (VU)T)  (6) 

In  addition  the  turbulent  stress  tensor  has  been  modeled  using  an  isotropic 
eddy  viscosity  such  that: 

uT  -  -  p  u'  u*  *  2uT  (7) 

where  Ut»  the  turbulent  viscosity,  is  determined  by  a  suitable  turbulence 
model.  Turbulence  modeling  is  described  in  some  detail  later. 

Equation  3  contains  a  mean  heat  flux  vector  defined  as  follows: 

q  -  -  K  (VT)  (8) 

and  a  turbulent  heat  flux  vector  defined  as: 


where  <  and  kt  are  the  molecular  (laminar)  and  turbulent  thermal 
conductivities,  respectively. 

Also  appearing  in  Eq.  3  is  the  mean  flow  dissipation  term  *. 

*  =  2v  0:  ID-  (|  p  -  Kj  (7  •  U)2  (10 


The  equation  of  state  for  a  perfect  gas 


p  =  PRT  (11 

where  R  is  the  gas  constant,  the  caloric  equation  of  state 

e  *  c  T  (12 

v 

and  the  definition  of  static  enthalpy 

h  =  c  T  (13 

P 

supplement  the  equations  of  motion. 


Finally  the  flow  properties  u,  <  and  Kg  are  determined  using  the  following 
constitutive  relations. 

The  molecular  viscosity  u  is  determined  using  Sutherland's  law. 
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o _ J_ 
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where  S|  =  100°K  for  air. 

The  bulk  viscosity  will  be  assumed  to  be  zero, 

Kg  -  0  (15 


and  the  thermal  conductivity  may  be  determined  by  use  of  a  relation  similar 
to  Sutherland's  law  viz. 
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where  So  =  194°K  for  air. 


Turbulence  Modeling 


Three  possible  turbulence  models  are  contained  within  the  disk  pump 
code.  These  are  ( i)  zero  equation  model  -  (mixing  length  model), 

(it)  one-equation  model  -  (a  turbulence  energy-algebraic  length  scale  model), 
and  ( l  i  i )  two-equation  model  -  (a  turbulence  energy-turbulence  dissipation 
model).  Among  these  models  (i)  and  (iii)  were  included  in  the  computation. 
Therefore,  a  detailed  explanation  would  be  given  to  only  (i)  and  (iii). 

Zero  Equation  Model  -  (Mixing  Length) 

Of  all  available  turbulence  models,  Prandtl's  mixing  length  model  is 
probably  still  the  most  widely  used.  The  model  was  originally  developed  for 
use  in  unseparated  boundary  layer  flow  situations  and  has  been  shown  to 
perform  well  under  such  conditions.  In  the  cases  described  in  this  report 
the  mixing  length  model  has  ben  used  during  the  early  transient  period  of 
flow  development  in  a  disk  pump.  An  advantage  of  the  method  from  the  point  of 
view  of  economy  is  that  is  does  not  require  additional  transport  equations  to 
model  the  effect  of  turbulence,  but  rather  relates  the  Reynold's  shear  stress 
to  mean  flow  quantities  via: 
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where 


where 
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where  d  is  the  normal  distance  to  the  nearest  wall  and  D  is  van  Driest 
damping  coefficient  given  by 

D  =  1  -  exp(-  y  /A  ) 


where  the  local  shear  stress  t £  is  obtained  from 


and 


is  defined  by  Eq.  6. 


t  =  (2  CD  :  ID  >1/2 
A 


One  problem  in  the  mixing  length  formulation  is  the  definition  of  6. 

In  boundary  layers  the  streamwise  velocity  u  approaches  an  edge  velocity, 
ue ,  asymptotically.  However,  a  monotonic  approach  to  an  asymptotic  edge 
velocity  is  not  characteristic  of  Navier-Stokes  solutions  particularly  in 
complicated  recirculating  flows  such  as  that  would  be  considered  in  the 
present  report.  In  order  to  avoid  the  problem  of  determining  the  boundary 
layer  edge,  6,  as  defined  in  the  usual  boundary  layer  context,  i.e.,  6  is  the 
distance  from  the  wall  at  which  u/ue  =  0.99,  the  following  relation  is 
used.  , 


2.0d 


(*l/ 9max~c) 


(20) 


In  other  words,  6,  is  taken  as  twice  the  distance  (measured  away  from  the 
nearest  wall)  for  which  q/qmax  -  c.  The  value  of  c  used  in  the  present 
effort  was  0.90. 


Two-Equation  Model  -  (k-e) 

The  k-e,  two-equation  turbulence  model  [6,  7]  in  which  both  the 
turbulence  kinetic  energy  and  the  turbulence  dissipation  rate  are  governed  by 
transport  equations  represents  a  more  general  model.  In  this  approach  the 
k-equation  is  as  given  by: 


3(pk)  +  V  .  (pUk)  =  V  •  Vk  +  2p  (  [D  :  CD  )  -  pe  -  2pv  (Vk1/2)2  (21) 
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and  the  e-equation  by: 
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where  following  Reference  6 


and 
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and  Rf  is  defined  as: 


(24) 


However,  attempts  to  solve  Eqs .  21  and  22  without  modification  fail 
because  an  appropriate  boundary  condition  for  e  at  a  solid  boundary  is 
difficult  to  prescribe  such  that  Eq.  22  is  satisfied.  In  order  to  circumvent 
this  problem,  Eq.  22,  the  turbulence  dissipation  equation,  has  been  modified 
by  the  inclusion  of  the  term: 

-  2uut  (V2U)2  (25) 

As  discussed  by  Jones  and  Launder  (Ref.  6)  inclusion  of  this  term  leads  to 
better  agreement  with  experimentally  measured  k  distributions  in  the 
near-wall  region.  In  addition,  Jones  and  Launder' s  (Ref.  6)  inclusion  of  the 
term : 

-  2pv  (Vk1/2)2  (26) 


in  Eq .  21  is  a  device  which  models  new  wall  dissipation  and  thus  allows  e=0 
to  be  prescribed  as  a  function  boundary  condition. 


Transformation  of  Governing  Equations 


The  governing  equations  for  the  present  problems  have  been  given  in  the 
previous  section  in  vector  form  (Eqs.  1-3  and  21  -  22).  However, 
implementation  of  a  solution  procedure  requires  that  these  equations  be 
transformed  into  an  appropriate  coordinate  system.  Therefore,  the  governing 
equations  written  in  a  cylindrical-polar  coordinate  system  are  transformed 
with  a  general  Jacobian  transformation  of  the  form 


yj  =  yj  (x^  t  x2  ,  X3  ,  t) 

T  =  t 


(27) 
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where  (xj,  X2,  X3)  -  (r,  0,  z)  are  the  original  cylindrical  polar 
coordinates.  The  velocity  components  remain  the  components  (Uj,  U2, 

U3)  in  the  (xj,  X2,  X3)  coordinate  directions,  respectively.  The  new 
independent  variables  yJ  are  the  computational  coordinates  in  the 
transformed  system. 

Application  of  the  Jacobian  transformation  requires  expansion  of  the 
temporal  and  spatial  derivatives  using  the  chain  rule,  i.e.. 


and 


where 


(28) 


(29) 


(30) 


The  relations  Eqs .  29  -  30  are  first  substituted  into  the  governing  equations 
written  in  cylindrical  polar  coordinates.  Then  the  resulting  equations  are 
multiplied  by  the  Jacobian  determinant  of  the  inverse  transformation. 
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and  the  equations  are  cast  into  a  conservative  form  using  the  following 
re  lat  ions 
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where 


yJ,i  s  JyJ»i 


(34) 


yJ  ,t  H  JyJ 


It  should  be  noted  that  Equation  33  expresses  a  geometric  conservation  law 
and  plays  an  important  role  in  enabling  the  equations  of  motion  to  be  cast 
into  conservative  form. 

The  particular  conservation  form  developed  implies  that  all  factors  involving 


the  radial  coordinate  r  *  x^  remain  as  they  were  before  the  Jacobian 
transformation.  The  resulting  equations  are  presented  below  as  Eqs.  37  -  46. 
The  geometric  relations  Eq.  32  -  33  may  be  obtained  from  the 
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transformation  relations  for  y,j|  and  y,^  in  terms  of  the  inverse 
transformat  ion  derivatives 
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rJ.t  =  -  I  yJ. 
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k  =  1 


The  transformed  governing  equations  may  be  written  in  the  following 


compact  form: 


3  .  3  . 

3(JW)  =  -  l  3  UyJ,,w)  -  l  ($.  3  (JyJ,.F  ) 

9T  J  =  1  3yJ  J  =  1  3yJ 


♦  Yi  J_  (Jyj*i  p£)  +  ?i  _2_  (Jy’.i  g.)) 


+  js  +  JC 


where 


y,J  =3jr 
t  3t 


y»Ji  H  hi 


Further,  the  coefficients  0{,  Yf,  Cj  are  given  by 


0  =  J_  ,  6  -  J. ,  8  -  1 

r  2  r  3 


Y  =  1,  Y  »  _1,  Y  -  1 

1  2  r  3 


1  ,  C  -  J.,  ?  =1 

m  2  r  3 


and  m  =  1  for  all  equations  except  the  X2  ~  direction  momentum  equation, 
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for  which  m  =  2.  The  vector  variables  used  in  Eq.  37  are  defined  as 
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where  n  =  1  for  i  *  1  and  n  =  0  for  i  =  2,  3. 


Note  that  the  velocity  components  (Oj,  U2,  U3)  are  the  cylindrical- 
polar  velocity  components,  and  t£j  is  the  stress  tensor  written  in 
cylindrical-polar  coordinates.  The  molecular  and  turbulent  stress  tensors 
may  be  written  as 


T . .  =  2p  ,,  D  . . 
ij  eff  ij 
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The  derivatives  required  in  Eqs .  43-44  must  be  expressed  in  terms  of  the 
computational  coordinates  yJ  using  the  chain  rule,  (Eq.  29). 

Finally,  the  vector  S  contains  source  terras  and  certain  differential 
terms  which  do  not  conform  to  the  basic  structure  of  Eq.  37,  and  the  vector 
contains  the  additional  curvature  terms  due  to  the  cylindrical-polar 
coordinate  system. 
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Numerical  Procedure 


The  numerical  procedure  used  to  solve  the  governing  equations  is  a 
consistently  split  linearized  block  implicit  (LBI)  scheme  originally 
developed  by  Briley  and  McDonald  (Ref.  4).  A  conceptually  similar  scheme  has 
been  developed  for  two-dimensional  MHD  problems  by  Lindemuth  and  Killeen 
(Ref.  9).  The  procedure  is  discussed  in  detail  in  Refs.  4  and  5.  The  method 
can  be  briefly  outlined  as  follows:  the  governing  equations  are  replaced  by 
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an  implicit  time  difference  approximation,  optionally  a  backward  difference 
or  Crank-Nicolson  scheme.  Terras  involving  nonlinearities  at  the  implicit 
time  level  are  linearized  by  Taylor  series  expansion  in  time  about  the 
solution  at  the  known  time  level,  and  spatial  difference  approximations  are 
introduced.  The  result  is  a  system  of  multidimensional  coupled  (but  linear) 
difference  equations  for  the  dependent  variables  at  the  unknown  or  implicit 
time  level.  To  solve  these  difference  equations,  the  Douglas-Gunn  (Ref.  10) 
procedure  for  generating  alternating-direction  implicit  (ADI)  schemes  as 
perturbations  of  fundamental  implicit  difference  schemes  is  introduced  in  its 
natural  extension  to  systems  of  partial  differential  equations. 

This  technique  leads  to  systems  of  coupled  linear  difference  equations  having 
narrow  block-banded  matrix  structures  which  can  be  solved  efficiently  by 
standard  block-elimination  methods. 

The  method  centers  around  the  use  of  a  formal  linearization  technique 
adapted  for  the  integration  of  initial-value  problems.  The  linearization 
technique,  which  requires  an  implicit  solution  procedure,  permits  the 
solution  of  coupled  nonlinear  equations  in  one  space  dimension  (to  the 
requisite  degree  of  accuracy)  by  a  one-step  noniterative  scheme.  Since  no 
iteration  is  required  to  compute  the  solution  for  a  single  time  step,  and 
since  only  moderate  effort  is  required  for  solution  of  the  implicit 
difference  equations,  the  method  is  computationally  efficient;  this 
efficiency  is  retained  for  multidimensional  problems  by  using  what  might  be 
termed  block  ADI  techniques.  The  method  is  also  economical  in  terms  of 
computer  storage,  in  its  present  form  requiring  only  two  time-level6  of 
storage  for  each  dependent  variable.  Furthermore,  the  block  ADI  technique 
reduces  multidimensional  problems  to  sequences  of  calculations  which  are 
one-dimensional  in  the  sense  that  easily-solved  narrow  block-banded  matrices 
associated  with  one-dimensional  rows  of  grid  points  are  produced.  A  more 
detailed  discussion  of  the  solution  procedure  as  discussed  by  Briley,  Buggeln 
and  McDonald  (Ref.  11)  is  given  in  the  Appendix  A. 


CURRENT  EFFORTS 


Objective 


A  primary  objective  of  the  current  efforts  was  to  demonstrate  the 
ability  of  a  modified  version  of  the  labyrinth  seal  code  to  be  applied  to  a 
disk  pumping  cavity  of  the  type  d  in  Fig.  1.  Since  the  present  Phase  1  work 
focused  upon  demonstration  and  assessment,  model  axisymmetric  cavity  shapes 
were  considered.  This  configuration  allows  swirling  flows  in  the  cavity 
although  the  computation  is  performed  in  two  dimensions  since  no  azimuthal 
variation  is  assumed.  The  current  type  of  disk  pump  operates  with  one  cavity 
wall  fixed,  while  the  other  rotating.  Air  enters  at  a  location  near  the 
inner  radius,  and  leaves  along  the  disk  periphery  as  shown  in  Fig.  2.  The 
disk  periphery  was  assumed  to  be  an  exit  boundary  in  the  sense  that  air  was 
not  permitted  to  enter  the  cavity  from  the  free  stream  and  similarly  flow 
within  the  cavity  was  not  permitted  to  leave  via  the  inlet  on  the  axis  of 
rotation.  As  shall  be  discussed,  the  specific  objectives  under  the  Phase  I 
effort  have  been  met  completely. 


Coordinate  System 


The  coordinate  system  is  an  important  component  of  the  Navier-Stokes 
analysis.  An  inappropriate  coordinate  system  may  lead  to  difficulty  in 
obtaining  a  converged  solution  and  even  exhibit  physically  unrealistic 
predictions  due  to  geometric  truncation  error.  Therefore,  generation  f  a 
viable  system  is  mandatory.  Any  coordinate  system  used  in  the  analysis 
should  satisfy  conditions  of  (i)  generality,  (ii)  smoothness, 

(iii)  resolvability  and  (iv)  allow  easy  application  of  boundary  conditions. 
Obviously,  a  coordinate  system  must  be  sufficiently  general  to  allow 
application  to  a  wide  class  of  problems  of  interest  if  it  is  to  be 
practical.  The  metric  data  associated  with  the  coordinte  system  must  be 
sufficiently  smooth  so  that  the  variation  from  grid  point  to  grid  point  does 
not  lead  to  a  numerical  solution  dominated  by  metric  coefficient  truncation 
error;  it  should  be  noted  that  this  requirement  differs  from  the  requirement 
of  the  existence  of  a  specified  number  of  transformation  derivatives. 

The  coordinate  system  must  resolve  flow  regions  where  rapid  flow  field 
change  occurs.  Finally,  coordinates  should  allow  accurate  implementation  of 
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boundary  conditions.  The  specification  of  boundary  surfaces  which  do  not 
fall  upon  coordinate  lines  or  at  specific  grid  points  may  present  a  difficult 
problem  for  Navier-Stokes  analyses.  The  problem  is  considerably  more  severe 
in  viscous  flows  where  no-slip  conditions  on  solid  walls  can  combine  with 
boundary  condition  truncation  error  to  produce  numerical  solutions  which  are 
both  qualitatively  and  quantitatively  in  error.  Thus,  the  property  required 
of  any  coordinate  system  to  be  used  in  a  viscous  analysis  should  be  that 
boundary  surfaces  fall  on  coordinate  lines.  Such  a  system  is  called  a 
body-fitted  coordinate  system. 

Many  alternative  approaches  are  available  for  generating  body-fitted 
coordinates,  although  in  general  they  can  be  categorized  into  three  generic 
groups.  These  are  conformal,  algebraic  and  elliptic.  Conformal  systems  are 
based  upon  well  known  conformal  transformation  techniques.  Algebraic  or 
constructive  techniques  are  typified  by  the  method  of  Eiseman  (Ref.  12)  while 
methods  using  elliptic  partial  differential  equations  are  in  the  main  based 
on  the  work  of  Thompson,  Thames  and  Mastin  (Ref.  13).  While  the  main 
advantage  of  the  use  of  elliptic  generation  is  that  the  method  is  applicable 
to  complex  geometries  without  the  need  to  modify  the  basic  approach  for  each 
new  conf igurat ion,  algebraic  methods  have  the  apparent  advantage  of 
conceptual  simplicity  for  complex  boundary  geometries.  In  the  present  work 
the  constructive  approach  has  been  adopted,  although  the  Thompson  approach  is 
presently  available  at  SRA  and  could  be  utilized  in  the  future  efforts 
associated  with  disk  pump  flow  analysis. 

In  the  case  of  plain  disk  cavity,  the  coordinate  system  is 
straightforward.  A  simple  cylindrical  coordiante  system  is  utilized  because 
of  its  axisymmetric  shape  as  shown  in  Fig.  3.  Horizontal  lines  represent  the 
axial  coordinate  lines,  while  the  radial  coordinate  lines  are  represented  by 
transverse  lines.  Partly  embedded  solid  regions  were  assumed  within  the 
domain  of  interest  for  coordinate  generation.  As  a  result,  reentrant  corners 
were  generated.  A  further  distribution  of  grid  points  within  the 
computational  domain  is  required  to  obtain  high  grid  resolution  in  regions 
where  rapid  flow  variations  are  expected.  This  is  accomplished  by  use  an 
analytical  transformation  technique  developed  by  Oh  (Ref.  14).  The 
transformation  function  is  composed  of  a  series  of  complementary  error 
functions.  One  of  the  advantages  in  the  use  of  this  transformation  function 
lies  in  the  capability  of  controlling  the  local  grid  points  distribution 
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without  effecting  the  distribution  at  other  locations  of  interest.  For 
further  details  of  the  technique,  the  reader  should  refer  to  Ref.  14. 

On  the  other  hand,  the  shaped  disk  cavity  grid  coordinate  system 
requires  a  calculation  or  specification  of  arc  lengths  on  the  boundaries. 
Then,  the  distribution  of  grid  points  on  and  within  the  boundaries  are 
obtained  based  on  the  parametrized  arc  length.  The  same  transformation 
technique  (Ref.  14)  was  used  to  pack  the  grid  points  in  the  region  of  large 
gradients  in  flow  variables.  The  final  coordinate  systems  used  in  the 
present  work  are  shown  in  Fig.  3. 

Navier-Stokes  Analysis 

Having  generated  a  viable  coordinate  system,  the  next  step  in  the 
process  is  to  obtain  a  solution  of  the  Navier-Stokes  equations  for  the 
specified  coordinate  system  and  flow  conditions.  The  analysis  was  performed 
in  the  same  ax isymmetr ic  configurations  as  shown  in  Fig.  2.  One  of  the 
cavity  walls  was  assumed  to  rotate  about  the  axis  at  3600  rpm,  while  the 
other  wall  was  stationary.  The  major  parameters  representing  the  turbine 
disk  test  conditions  are  axial  and  tangential  Reynolds  numbers.  The  axial 
Reynolds  number  is  defined  as  Rez  =  ra/irur0  where  m  is  the  mass  flow  rate 
through  the  inlet,  u  is  the  laminar  viscosity  and  rQ  is  the  maximum  radius 

•y 

of  rig,  while  the  tangential  Reynolds  number  is  defined  as  Ret  =  p0ro 

U 

there  p  is  the  density,  0  is  the  angular  velocity  of  rotating  disk  wall.  The 
axial  and  tangential  Reynolds  number  selected  in  the  present  analysis  was 
2.2  x  104,  and  4.0  x  10^,  respectively.  The  laminar  viscosity,  p,  in  the 
present  analysis  was  1.905  x  10“  (Newton-sec/m  ) .  Other  flow  conditions 
were  derived  from  both  Reynolds  numbers  and  laminar  viscosity.  Reference 
Mach  number  was  M  =  0.  135,  while  the  static  pressure  was  estimated  as 

2.3197  x  10  (Newton/m  ).  The  static  temperature  of  the  stream  was  312. 4*K. 

Meanwhile,  the  static  pressure  at  the  exit  boundary  was  assumed  to  be 
atmospheric.  Thus,  the  ratio  of  exit  static  pressure  to  the  inlet  pressure 
was  2.  Reference  length  used  in  the  computation  was  the  maximum  radius  of 

the  rig  which  was  0.2794m.  The  same  flow  conditions  were  used  for  both  plain 

and  shaped  disk  pump  cavity.  To  analyze  the  flow  in  both  cavities  101  mesh 
points  were  distributed  in  the  radial  direction.  In  the  axial  direction,  60 
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mesh  points  were  used  for  the  plain  cavity  case,  and  71  mesh  points  were  used 
for  the  shaped  cavity  case.  The  grid  points  were  tightly  packed  in  the 
vicinity  of  solid  walls  to  achieve  the  maximum  resolution  in  the  boundary 
layers.  For  both  cases  the  continuity,  three  momenta  equations  and  the 
constant  stagnation  enthalpy  condition  were  used  as  governing  equations. 
However,  because  of  the  axisymmetric  flow  configuration,  (no  change  of  flow 
variables  in  the  azimuthal  direction),  the  computation  required  only  two 
directions  (axial  and  radial  direction). 


Boundary  Conditions 

A  major  factor  in  obtaining  efficient  solutions  for  the  Navier-stokes 
equations  is  specification  of  appropriate  boundary  conditions.  Specification 
of  boundary  conditions  at  wall  is  relatively  straightforward,  but  proper 
specification  at  inflow  and  outflow  boundaries  presents  a  more  difficult 
problem.  The  present  effort  follows  the  work  of  Briley  and 

McDonald  (Ref.  15)  who  suggest  examining  the  charactersistics  of  the  inviscid 
problem  for  guidance  at  inflow  and  outflow  boundaries.  Appropriate  boundary 
conditions  which  were  set  in  the  present  calculations  are  as  follows: 


(i)  Inlet  (subsonic  inflow)  boundary  - 

Function  condition  on  uj-velocity  (radial  velocity) 

Two-layer  model  on  U3~velocity  (axial  velocity)  and  static  enthalpy 
Second  derivative  of  static  pressure  set  to  zero. 

(li)  Exit  (subsonic  outflow)  boundary  - 

Second  derivative  of  ui,  U3  set  to  zero 
Static  pressure  specified. 

(iii)  Wall  boundary  - 

No-slip  on  ui  and  U3 

Normal  (to  the  wall)  momentum  equation 

Function  condition  on  u  (zero  on  stationary  wall  and  nonzero  on 
rotat ing  wall) . 


.*  Vj  J*  Ip  _  j  ^  ± 


*  .'3*  PJ>*.  V 
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The  so-called  two-layer  model  used  at  the  inflow  boundary  is 
essentially  a  total  pressure  boundary  condition  applied  to  the  core  flow  with 
a  specified  boundary  layer  profile  shape  for  the  wall  region.  Matching  the 
static  pressure  at  the  edge,  defined  by  the  first  computational  point  from 
the  wall  at  which  1 113 1 /  | U3 1 max  was  greater  than  or  equal  to  0.99  on  the 
previous  time  step,  enables  calcuation  of  |u3|  at  this  point.  This  provides 
the  required  normalizing  value  for  the  pre-spec  if ied  boundary  layer  profile 
shape.  Overall,  the  method  provides  a  mechanism  for  drawing  mass  into  the 
system  in  order  to  satisfy  the  downstream  pressure  given  an  upstream  core 
total  pressure.  This  specification  corresponds  to  the  usual  wind 
tunnel  experiment  where  stagnation  conditions  are  set  in  an  upstream 
reservoir  and  static  pressure  is  set  at  some  downstream  location. 
Specification  of  an  inlet  mass  flux  could  be  accomplished  indirectly  by 
varying  the  downstream  static  pressure. 

Calculation 

Effective  running  of  the  disk  pumping  cavity  problem  requires  a  case 
running  protocol  to  obtain  rapid  convergence.  The  basic  strategy  for  running 
the  code  consists  of  several  steps.  First,  initial  conditions  within  the 
both  cavities  assumed  stagnant  flow.  Then,  as  the  specified  exit  pressure 
at  the  outflow  boundary  is  lowered,  mass  flow  is  drawn  through  the  inlet. 
Swirl  was  assumed  to  be  zero,  and  a  mixing  length  turbulence  model  was  used 
initially  to  develop  the  basic  flow  pattern  within  the  cavities  without  the 
complexities  of  a  swirl  equation  or  a  two-equation  turbulence  model. 

Then  two  different  methods  were  implemented  for  the  plain  and  shaped  disks. 

In  the  case  of  the  plain  disk,  rotation  effects  was  introduced  before  the 
two-equation  turbulence  model  replaced  the  mixing- length  model.  On  the  other 
hand,  in  the  case  of  the  shaped  disk,  the  order  was  reversed,  i.e.  the 
two-equation  turbulence  model  was  introduced  before  rotation  effects. 

Both  run  protocols  were  successful  and  neither  appeared  to  have  a  definite 
advantage  over  the  other.  In  both  cases  two-equation  turbulence  model 
equations  were  initially  solved  with  invariant  mean  flows  to  develop 
approximate  turbulence  energy  and  dissipation  fields  prior  to  coupling  these 
fields  with  the  mean  flow.  Then,  finally  all  equations  including  the 
turbulence  equations  were  solved  simultaneously  until  a  converged  solution 
was  obtained. 


Although  the  present  Navier-Stokes  analysis  can  be  used  both  for 
time-dependent  and  steady-state  flow  situations,  the  focus  of  the  present 
investigation  is  the  steady  state  solution.  In  such  a  case,  it  is  not 
necessary  to  accurately  follow  the  transient  motion  and  indeed  it  may  be 
advantageous  from  a  computational  efficiency  view  point  not  to  follow  the 
transient  motion  accurately,  if  this  accelerates  convergence  to  the  steady 
state.  The  present  approach  utilizes  the  matrix  conditioning  technique  of 
Refs.  16  and  17  to  accelerate  convergence  to  a  steady  state.  Using  the 
techniques  described  in  Refs.  16  and  17,  the  convergence  took  250 
time  steps  for  both  plain  and  shaped  disk,  respectively.  The  CRAY  CPU  time 
required  for  the  convergence  of  both  cases  was  less  than  2100  secs. 

Results 

The  calcuated  velocity,  Mach  number  and  stream  function  contours  of  two 
test  cases  are  presented  in  Fig.  4  through  Fig.  10.  Since  experimental  data 
is  not  available  for  direct  comparison,  discussion  will  primarily  focus  on 
qualitative  description  of  the  essential  physical  features  predicted  by  the 
computation.  Although  both  dependent  variable  contours  and  flow  streamlines 
are  presented,  the  clearest  picture  of  the  generated  flow  can  be  obtained 
from  the  flow  streamlines. 

In  the  case  of  a  plain  disk  pump,  a  large  recirculating  fluid  motion  is 
generated  within  the  cavity  as  shown  in  Figs.  8  and  9.  The  flow  pattern  in 
the  absence  of  wall  rotation  is  basically  determined  by  the  convection  of  the 
fluid  due  to  the  pressure  difference  between  inlet  and  exit.  In  the  present 
work,  the  static  pressure  at  the  exit  is  one-half  of  the  pressure  at  the 
inlet.  When  the  effects  of  disk  wall  rotation  are  considered,  the  flow 
pattern  change  is  apparent  including  the  radial  displacment  of  the 
recirculation  region  within  the  cavity  as  seen  in  Fig.  9.  However, 
considering  that  the  basic  flow  structure  still  remains  the  same  with  and 
without  rotation,  it  is  likely  that  centrifugal  pumping  action  is  not 
dominant  over  the  effects  of  convection  because  of  a  relatively  large 
pressure  drop  in  this  case. 

The  flow  pattern  within  the  shaped  disk  cavity  in  the  absence  of  wall 
rotation  is  different  from  the  plain  case.  It  is  likely  that  this  difference 
is  strongly  dependent  on  the  geometry  (Fig.  9-10).  Under  the  same 


nonrot.it  ing  and  pressure  condition  as  the  plain  disk  case,  the  flow  within 
the  shaped  disk  contains  only  a  small  recirculation  region  in  the  vicinity  of 
reentrant  corner  near  the  inlet.  When  the  effects  of  rotation  are 
considered,  the  flow  within  the  cavity  changes  substantially  as  shown  in 
Fig.  10.  The  centrifugal  pumping  action  causes  a  considerable  change  in  the 
flow  fieLd  for  the  shaped  disk  pump  case.  Since  such  a  complex  fluid  motion 
can  be  associated  with  many  physical  mechanisms,  it  may  be  too  premature  to 
draw  a  definite  conclusion  before  further  detailed  parametric  studies  are 
performed.  However,  it  is  likely  that  centrifugal  pumping  action  sweeps  the 
fluid  away  in  the  radial  outward  direction  near  the  rotating  cavity  wall  and 
mduces  a  large  recirculating  fluid  motion.  Velocity  and  Mach  number 
contours  as  shown  in  Figs.  4-7  present  the  further  details  of  flow  structures 
associated  with  both  disk  pump  shapes  for  the  cases  which  include  rotation. 

Of  particular  interest  is  the  contours  of  swirling  velocity  component  (V). 

The  figure  indicates  that  the  convection  is  closely  related  to  the  faster 
transfer  of  the  swirling  component  of  momentum  in  the  axial  direction  for 
both  shapes  of  disk. 

A  major  application  of  the  rotor  disk  is  to  pump  the  cooling  air  out 
into  the  hot  gas  path.  However,  since  the  disk  wall  (rotating)  is  exposed  to 
hot  temperature  environment,  the  fluids  within  the  cavity  especially  in  the 
vicinity  of  the  rotating  disk  may  heat  up  due  to  the  heat  transfer  through 
the  disk  wall.  As  a  result,  not  only  cooling  efficiency  loss  but  also 
structural  damage  may  arise.  Comparison  of  flow  pattern  in  both  shapes 
indicate  that  the  shaped  disk  is  more  efficient  in  preventing  such  a  hot 
temperature  region  developing  adjacent  to  the  disk  wall.  In  this  regard  it 
should  be  noted  that  the  closeness  of  the  streamlines  indicates  the  speed  of 
the  flow.  The  contours  in  Fig.  8  clearly  show  a  large  recirculation  region 
along  the  right  side  wall  for  the  plain  disk.  Fluid  trapped  in  this  region 
will  heat  up  and  lose  its  cooling  effectiveness.  Hence,  the  right  hand  wall 
in  the  plain  disk  would  be  expected  to  be  a  region  of  high  wall  temperature. 
In  the  case  of  the  rotating  shaped  disk,  the  situation  is  much  more  favorable 
as  shown  in  Fig.  8.  Here  cooling  problems  may  occur  on  the  upper  wall, 
however,  both  the  right  and  left  side  walls  appear  to  be  in  contact  with 
relatively  rapidly  moving  cooling  fluid.  It  should  be  noted  that  rotation 
makes  a  major  difference  in  the  flow  and  heat  transfer  patterns  for  the 
shaped  cavity.  In  view  of  the  strong  dependence  of  flow  field  within  the 
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cavity  on  the  shape  of  the  disk  wall,  various  geometric  shapes  should  be 
explored  to  find  an  optimum  shape  of  the  disk  pump. 

ESTIMATES  OF  TECHNICAL  FEASIBILITY 

Under  the  present  effort,  an  existing  Navier-Stokes  solver  developed 
for  the  analysis  of  labyrinth  seal  has  been  extended  and  applied  to 
demonstrate  its  capability  for  the  disk  pumping  problem.  Coordinate  systems 
were  generated  by  the  constructive  approach,  and  used  successfully  for  both 
plain  and  shaped  disk  pumps.  However,  a  more  general  elliptic  solver  has 
been  developed  at  SRA  to  generate  the  coordinate  systems  for  the  complex 
geometry  of  the  advanced  seal.  This  coordinate  generation  code  would  be 
available  for  the  future  analysis  of  disk  pump  problems  under  the  more 
complex  geometric  configuration  (Fig.  11). 

The  Navier-Stokes  code  was  used  to  calculate  the  disk  pump  flow 
fields  for  two  cases.  In  order  to  reduce  run  time  vectorization  of  the  code 
was  initiated  under  the  current  effort.  The  partially  vectorized  code 
reduced  run  times  by  factor  of  2  when  continuity  equation  and  two  momentum 
equations  were  solved.  Even  with  an  additional  momentum  equation,  the  run 
time  was  saved  by  25%  on  CRAY-1  computer.  Based  upon  other  efforts  at  SRA, 
it  is  expected  that  full  vectorization  would  reduce  run  times  by  a  further 
factor  of  5  to  10.  Although  in  the  present  work  a  constant  s'«»gration 
enthalpy  relation  was  assumed,  the  current  code  does  include  an  energy 
equation  (although  not  optimized)  and  has  been  successfully  used  in  heat 
transfer  calculations  (Ref.  18).  With  a  partially  vectorized  code,  the  work 
done  under  the  present  effort  has  allowed  steady  disk  pump  flow  fields  to  be 
generated  in  less  than  2100  secs  of  CRAY  CPU  time,  i.e.,  approximately 
$1200  at  a  commercially  available  rate. 

The  present  effort  has  concentrated  on  the  analysis  of  flow  fields 
within  the  modeled  cavity  shapes.  Because  of  the  general  nature  of  the 
Navier-Stokes  code  more  complex  disk  pump  configurations  can  be  considered 
in  the  future  (Fig.  11).  For  these  cases  the  elliptic  coordinate  system 
solver  would  be  used  to  generate  the  more  complex  configurations.  Results 
obtained  for  flow  within  both  cavities  are  very  encouraging  in  verifying  the 
capability  of  the  Navier-Stokes  procedure  for  such  a  practically  important 
problem.  Although  not  pursued  in  the  present  effort,  the  code  cold  be  used 
to  predict  the  heat  transfer  through  the  disk  pump  wall  by  solving  energy 
equat ion. 
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The  calculations  completed  under  the  Phase  I  effort  gave  physically  realistic 
flow  patterns,  and  showed  the  effect  of  cavity  shape  and  wall  rotation  upon 
the  flow  field.  The  run  times  used  were  encouraging,  and  with  full 
vector izat ion  it  is  expected  that  this  code  could  be  used  as  a  major  tool  for 
the  understanding  and  in  the  design  of  disk  cavity  configurations  including 
wall  heat  transfer  effects. 
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Fig.  8  -  Stream  Function 


Sketch  of  Compressor  Drum  Model 


APPENDIX  A  -  Split  LBI  Algorithm 


Linearization  and  Time  Differencing 

The  system  of  governing  equations  to  be  solved  consists  of  three/four 
equations:  continuity  and  two/ three  components  of  momentum  equation  in 
three/four  dependent  variables:  p,  u,  v,  w.  This  system  of  equations  can  be 
written  at  a  single  grid  point  in  the  following  form: 

3H($)/3t  *  D(*)  +  S(<|>)  (A-l) 

where  $  is  the  column-vector  of  dependent  variables,  H  and  S  are  column- 
vector  algebraic  functions  of  $,  and  D  is  a  column  vector  whose  elements  are 
the  spatial  differential  operators  which  generate  all  spatial  derivatives 
appearing  in  the  governing  equation  associated  with  that  element. 

The  solution  procedure  is  based  on  the  following  two-level  implicit 
time-difference  approximations  of  (Eq.  1): 

(Hn+1-  Hn)/At  *  0(Dn+1  +  S"+1)  +  (1-3)  (D°  +  S")  (A-2) 

where,  for  example,  Hn+1  denotes  and  At  *  tn+*  -  tn.  The 

parameter  3  (0.5  -  3  -  1)  permits  a  variable  time-centering  of  the  scheme, 
with  a  truncation  error  of  order  [At^,  (3  -  1/2)  At]. 

A  local  time  linearization  (Taylor  expansion  about  $n)  of  requisite 
formal  accuracy  is  introduced,  and  this  serves  to  define  a  linear 
differential  operator  L  such  that 


(A-3) 


x 


Similarly, 


Hn+1  =  H°+  (3H/3$)n  (<|in+1  -  <|>n)  +  0  (At2) 


(A-4) 


Sn+1  =  Sn+  (3S/3<J>)n  (<j>n+1  -  4»n)  +  0  (At2)  (A-5 ) 

Eqs.  (3-5)  are  inserted  into  Eq.  (2)  to  obtain  the  following  system  which  is 
linear  in  <frn+l 

(A  -  0At  Ln)  (<|>n+1  -  4>n)  =  At  (Dn  +  Sn)  (A-6) 

and  which  is  termed  a  linearized  block  implilcit  (LBI)  scheme.  Here,  A 
denotes  a  matrix  defined  by 

A  =  (3H/3<fr)n  -  0At  ( 3S/ 3$)n  (A-7) 

Eq.  (6)  has  0  (At)  accuracy  unless  H  =  $,  in  which  case  the  accuracy  is  the 
same  as  Eq.  (2). 


Special  Treatment  of  Diffusive  Terms 

The  time  differencing  of  diffusive  terms  is  modified  to  accomodate 
cross-derivative  terms  and  also  turbulent  viscosity  and  artificial  dissipa¬ 
tion  coefficients  which  depend  on  the  solution  variables.  Although  formal 
linearization  of  the  convection  and  pressure  gradient  terms  and  the  resulting 
implicit  coupling  of  variables  is  critical  to  the  stability  and  rapid  con¬ 
vergence  of  the  algorithm,  this  does  not  appear  to  be  important  for  the 
turbulent  viscosity  and  artificial  dissipation  coefficients.  Since  the 
relationship  between  ue  and  dj  and  the  mean  flow  variables  is  not  conven¬ 
iently  linearized,  these  diffusive  coefficients  are  evaluated  explicitly  at 
tn  during  each  time  step.  Notationally,  this  is  equivalent  to  neglecting 
terms  proportional  to  3  Pe/3<J>  or  3d j / 3 ♦  in  Ln,  which  are  formally  pre¬ 
sent  in  the  Taylor  series  expansion  (3),  but  retaining  all  terms  proportional 
to  ue  or  dj  in  both  Ln  and  Dn. 
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It  has  been  found  through  extensive  experience  that  this  has  little  if  any 
effect  on  the  performance  of  the  algorithm.  This  treatment  also  has  the  added 
benefit  that  the  turbulence  model  equations  can  be  decoupled  from  the  system  of 
mean  flow  equations  by  an  appropriate  matrix  partitioning,  and  solved  separately 
in  each  step  of  the  ADI  solution  procedure.  This  reduces  the  block  size  of  the 
block  tridiagonal  systems  which  must  be  solved  in  each  step  and  thus  reduces  the 
computational  labor. 

In  addition,  the  viscous  terms  in  the  present  formulation  include  a  number 
of  cross-derivative  terms  implicitly  within  the  ADI  treatment  which  follows,  it 
is  not  at  all  convenient  to  do  so;  and  consequently,  all  cross-derivative  terms 
are  evaluated  explicitly  at  tn.  For  a  scalar  model  equation  representing 
combined  convection  and  diffusion,  it  has  been  shown  by  Beam  and  Warming  that  the 
explicit  treatment  of  cross-derivative  terms  does  not  degrade  the  unconditional 
stability  of  the  present  algorithm.  To  preserve  notational  simplicity,  it  is 
understood  that  all  cross-derivative  terms  appearing  in  Ln  are  neglected  but 
are  retained  in  Dn.  It  is  important  to  note  that  neglecting  terms  in  Ln  has 
no  effect  on  steady  solutions  of  Eq.  (6),  since  5  0  and  thus  Eq.  (6) 

reduces  to  the  steady  form  of  the  equations:  Dn  +  Sn  =  0.  Aside  from 
stability  considerations,  the  only  effect  of  neglecting  terms  in  Ln  is  to 
introduce  an  0  (At)  truncation  error. 

Consistent  Splitting  of  the  LBI  Scheme 

To  obtain  an  efficient  algorithm,  the  linearized  system  (A-6)  is  split  using 
ADI  techniques.  To  obtain  the  split  scheme,  the  multidimensional  operator  L  is 
rewritten  as  the  sum  of  three  "one-dimensional"  sub-operators  Lj_  (i  =  1,  2,  3) 
each  of  which  contains  ail  terms  having  derivatives  with  respect  to  the  i-th 
coordinate.  The  split  form  of  Eq.  (6)  can  be  derived  either  as  in 
(Refs.  4  and  5)  by  following  the  procedure  described  by  Douglas  and  Gunn 
(Ref.  10)  in  their  generalization  and  unification  of  scalar  ADI  schemes,  or  using 
approximate  factorization.  For  the  present  system  of  equations,  the  split 
algorithm  is  given  by 


(A  -  BAtLj)  ($ 
(A  -  BAtL^)  (♦ 
(A  -  BAtLj)  ($ 


*  -  $")  -  At  (Dn  +  Sn)  (A-8) 

**  -  ♦“)  -  A  (♦* **  -  ♦")  (A-9) 

0+1  -  ♦")  -  A  (***  -  ♦")  (A- 10) 


where  and  $**  are  consistent  intermediate  solutions.  If  spatial 
derivatives  appearing  in  and  D  are  replace  by  three-point  difference 
formulas,  as  indicated  previously,  then  each  step  in  Eqs.  (8-10)  can  be  solved  by 
a  block-tridiagonal  elimination. 
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